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LITERATURE  SURVEY 


THE  ORIGINS  OF  LUBRICATION  THEORY 

After  Towers’  (1)  experimental  findings  concerning  the  existence  of 
pressures  significantly  above  atmospheric  pressure  in  the  oil  film  of  a 
journal  bearing,  Reynolds  (2)  developed  equations  to  predict  the  pressure 
distribution  due  to  the  viscose  forces.  This  was  the  first  theoretical 
explanation  of  a  practical  phenomena  which  hud  been  unknowingly  used  for 
years . 

To  east?  the  problem  of  boundary  conditions  required  for  the  solution  of 
Reynolds  equations,  Ocvirk  (3)  produced  a  short  bearing  approximation  with 
a  180  degree  positive  pressure  arc  and  circumferential  flow  neglected  and 
Sommer fe Id  produced  a  long  bearing  approximation  by  assuming  negligible 
pressure  variation  in  the  axial  direction.  These  works  formed  the 
foundation  for  later  developments  in  lubrication  theory. 


EARLY  DEVELOPMENTS  IN  ROTOR  DYNAMICS 

In  1869  Rankine  (4)  mathematically  examined  the  equilibrium  conditions  of  a 
mass  mounted  on  a  flexible  rotating  shaft  supported  in  frictionless 
bearings.  Due  to  his  neglect  of  Coriolis’  forces,  he  concluded  that  the 
rotor  is  always  unstable  above  the  shaft  critical  speed.  Later  in  1885  de 
Laval  experimentally  demonstrated  that  super  critical  operation  was 
possible  but  manufacturers  experienced  many  design  failures  due  to 
excessive  vibration  at  supercritical  speeds.  Dunkerley  (7)  in  1894  using  a 
multi-mass,  multi-stiffness  model  of  a  shaft  system,  theoretically 
predicted  the  ranges  of  "critical  speeds"  characterised  by  rough  running 
and  large  force  transmission  to  the  foundations. 
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Jeffcott  (5)  justified  the  use  of  an  elegantly  simple  single  mass,  flexible 
rotor  model  by  revealing  the  mechanism  of  shaft  resonance.  He  concluded 
that  below  the  shaft  pin-pin  critical  speed  the  rotor  centre  precesses  at 
the  same  angular  velocity  and  in  phase  with  the  shaft  imbalance.  This  was 
termed  "syncronous  precession".  At  the  shaft  critical  speed  this  phase 
difference  was  shown  to  be  90*  and  well  above  the  critical  speed  the 
unbalance  was  180* out  of  phase  with  the  shaft  centre.  The  shaft  was  shown 
to  reduce  its  flexure  amplitude  above  critical  speed  due  to  this 
"inversion"  process  (shaft  centre  line  rotating  about  its  mass  centre) 
which  led  in  the  1920 ’s  to  the  use  of  supercritical  rotors,  termed 
flexible. 

However,  the  problem  of  low  frequency  shaft  whipping  whilst  operating  at 
supercritical  speeds  remained  a  problem  and  so  in  1924  Newkirk  (6)  was 
commissioned  to  experimentally  investigate  the  problem.  Using  both 
built-up  and  non-built-up  rotors  supported  on  rolling  element  bearings  he 
noted  that  the  onset  speed  of  whirling  and  the  whirl  amplitude  were 
unaffected  by  rotor  unbalance  and  always  exhibited  the  same  precession 
speed.  The  threshold  speed  could  be  increased  by  increasing  the  foundation 
flexibility,  introducing  foundation  damping,  misaligning  the  bearings  and 
by  applying  axial  thrust  to  the  bearings.  The  whirling  was  found  only  to 
occur  with  built-up  rotors.  By  introducing  internal  friction  into 
Jeffcotts’  model,  Newkirk  showed  that  the  attenuating  effect  of  internal 
friction  became  an  exciting  force  well  above  critical  speed. 

I’oritsky  (8)  investigated  the  phenomena  using  a  similar  flexible  rotor 
model  but  contrary  to  Newkirk’s  findings  he  concluded  that  the  introduction 
of  system  foundation  flexibility  lowers  the  critical  speed  of  the  system 
and  hence  reduces  the  stability  threshold. 
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The  importance  of  foundation  flexibility  and  damping  (impedance)  on  the 
systems  dynamic  behaviour  had  by  this  time  been  realised  but  the 

theoretical  and  experimental  contradictions  as  to  the  effect  of  changing 
the  support  impedance  or  the  reasons  for  such  effects  were  unresolved. 

Work  by  Stodola  (10)  in  1925  showed  that  the  presence  of  the  oil  film 

journal  bearings  in  an  otherwise  rigid  system  would  reduce  the  system’s 
critical  speed  to  well  below  that  calculated  for  the  rigidly  supported 
system.  As  a  consequence  of  this  many  authors  (11,  12,  14-18)  have 

attempted  to  accurately  model  the  oil  film  effect  on  system  response  by 
developing  stiffness  and  damping  coefficients  to  represent  the  oil  film  in 
the  journal  bearing. 

Gunter  (13)  in  1966  further  developed  Jeffcotts’  model  of  a  single  point 
mass  flexible  rotor  with  bearing  radial  stiffness  and  internal  friction  to 
show  the  effect  of  foundation  flexibility  effects.  The  interaction  of 
support  stiffness  assymetry,  rotor  stiffness  and  damping  factors  is 
considered.  From  this  analysis  he  showed  that  for  a  given  value  of 
external  damping  there  is  a  value  of  support  flexibility  which  will  make 
the  system  entirely  stable  for  all  speeds.  In  the  absence  of  damping, 
assymetric  support  stiffness  is  shown  to  improve  stability  above  that  of 
the  symmetric  support.  However  this  results  in  a  non-circular  procession 
loci  of  the  shaft  centre  in  the  bearing  plane  with  the  associated  increase 
in  transmitted  force  at  the  ellipse  tails.  The  symmetric  support  produces 
circular  loci  (in  the  stable  syncronous  zone)  and  with  the  addition  of 
correct  support  damping,  increases  the  stability  threshold  to  its  highest 
value. 

The  assumption  of  negligible  gyroscopic  forces  is  generally  accepted  but 
the  single  mass  model  does  limit  the  models  representation  of  practical 
rotor  systems.  No  new  experimental  results  are  provided  but  Gunter  does 


compare  his  theoretical  results  with  the  experimental  work  of  Newkirk  and 
of  Pinkus  (7)  and  explains  some  of  the  apparent  anomalies  of  their 
findings. 


DOUBLE  BEARINGS 

Some  on  the  first  work  on  "double  bearings"  was  carried  out  by  Shaw  and 
Nussdorfer  in  1947  (19).  They  developed  a  form  of  double  journal  bearing 
with  a  floating  ring  in  between  the  two  oil  films.  The  object  of  the 
bearing  was  to  reduce  the  high  losses  associated  with  the  usual  single  film 
bearing  and  to  increase  oil  flow  to  allow  cooler  running  conditions.  The 
shaft  was  centralised  by  a  cantilever  mechanism  and  the  floating  centre 
ring  was  free  to  rotate.  Unfortunately,  despite  achieving  their  design 
objectives,  they  encountered  severe  instability  of  the  bearing  and  startup 
ring  seizure  when  the  bearing  was  statically  loaded,  a  problem  which 
remained  unsolved. 

In  196H  Orcutt  and  Ng  (20)  investigated  a  double  (floating  ring)  bearing 
similar  to  that  of  Shaw  and  Nussdorfer.  In  particular  they  investigated 
the  effect  of  supply  pressure,  relative  film  clearance  ratio  and  lubricant 
viscosity  on  bearing  performance.  Whirl  instability  within  the  design 
operating  range  again  limited  the  number  of  results  taken.  However 
stability  was  improved  by  bearing  oil  flow  restriction,  viscosity  increase, 
avoidance  of  excessively  high  supply  pressures  and  by  lowering  the  inner  to 
outer  film  clearance  ratio  from  1.3  to  0.7  (shown  by  comparing  whirl 
ampl i tudes) . 

The  stability  improvements  were  not  however  of  sufficient  magnitude  to 
completely  solve  the  instability  problems  encountered. 


5. 


SQUEEZE  FILM  DEVELOPMENT 


To  vary  the  support  characteristics  of  rotor  systems,  squeeze  films  can  be 
introduced  around  the  hearing  system  used.  The  stiffness  and  damping 
properties  of  the  oil  film  are  used  to  attenuate  the  vibrational  amplitudes 
and  the  forces  transmitted  to  the  foundations,  and  to  increase  the 
stability  threshold  speed  of  the  system.  However  these  factors  are 
sometimes  conflicting. 

There  are  two  main  design  techniques  used.  First,  t  he  spring  centralised 
squeeze  film  which  uses  a  squirrel  cage  type  spring  or  similar  to  overcome 
the  shaft  gravity  load  and  hence  provide  zero  eccentricity  in  the  squeeze 
film  at  zero  shaft  speed;  the  rotor  consequently  acts  as  a  vertical  rotor. 
The  second  type  is  the  non-central ised  squeeze  film  which  relies  on  shaft 
imbalance  which,  when  this  force  exceeds  the  gravity  load,  lifts  the  main 
bearing  in  the  squeeze  film.  Subsequently  the  dynamic  stiffness  of  the 
squeeze  film  results  in  the  lift  being  maintained.  This  second  form  tends 
to  be  cheaper  and  useful  when  space  is  restricted.  In  both  types  the  inner 
squeeze  film  surface  is  prevented  from  rotating  by  dogs  or  a  similar  device 
but  the  inner  surface  is  free  to  move  within  the  oil  film  in  order  to 
generate  the  film  forces  required. 

In  1979  Nilolajsen  and  Holmes  (30)  experimentally  examined  a  squeeze  film 
sup  orting  a  hydrodynamic  bearing  and  shaft.  The  separating  ring  was 
prevented  from  rotating  and  supported  by  three  Belleville  washer  stacks  so 
that  the  film  could  be  centralised  and  the  support  stiffness  could  be 
varied  independently  to  the  support  damping  which  was  changed  by  the  use  of 
various  viscosity  oils  and  squeeze  film  clearances. 
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Significant  reductions  in  response  amplitudes  were  obtained  by  increasing 
support  damping  but  the  critical  speed  of  the  system  showed  little 
variation  with  support,  stiffness  in  experiments  despite  theoretical 
predictions  to  the  contrary.  This  work  shows  how  attenuation  of  amplitude 
of  vibration  can  be  achieved  by  selecting  the  correct  value  of  support 
damping  as  discussed  in  Gunter’s  work  (13). 

The  dynamics  of  a  simple  single  mass  rotor  rigid  shaft  with  squeeze  film 
supported  rolling  element  bearings  was  analysed  using  Ocvirks’  solution  of 
Heynolds  equation,  by  Cooper  (21).  lie  predicted  three  orbit  solutions  for 
a  given  unbalance.  White  (22)  experimentally  verified  the  predicted  orbits 
showing  the  higher  eccentricity  (whirl)  mode  and  the  lower  eccentricity 
( inverted)  mode  which  were  stable  and  the  intermediate  solution  which  has 
always  been  unstable.  The  "jump"  speeds  (speeds  at  which  the  bearing  tends 
to  jump  from  the  lower  eccentricity  mode  to  the  upper  eccentricity  mode) 
were  shown  by  White  to  be  as  predicted  by  Cooper.  However  from  the 
correlation  of  White’s  experimental  and  theoretical  results  it  is  shown 
that  White’s  model  is  only  accurate  for  eccentricity  ratio  less  thun  0.4. 

SQUEEZE  FILM  SUPPOKTEU  HO LLING  ELEMENT  BEARINGS 

In  the  1970’s  work  was  focused  on  the  prevention  of  large:  whirl  amplitudes 
exhibited  by  rotating  shafts  mounted  in  rolling  element  bearings  by 
providing  additional  support  damping  in  the  form  of  squeeze  film 
surrounding  the  bearing.  The  excitation  of  this  system  by  the  rolling 
elements  has  been  shown  by  several  authors  (26,  47)  to  be  negligible. 

In  1974  Mohan  and  Hahn  (23)  theoretically  examined  the  centrally  preloaded 
squeeze  film  damper  system  for  rigid  rotors  to  predict  the  squeeze  film 
forces  generated  and  their  effect  on  system  response.  Using  the  short 
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bearing  approximation  of  Reynolds  equation,  the  film  forces  and  hence  the 
velocity  (damping)  forces  and  deflection  (stiffness)  forces  were  used  to 
define  the  linear  stiffness  and  damping  coefficients  of  the  film.  The 
model  showed  the  jump  from  uninverted  to  inverted  mode,  predicted  by  Cooper 
and  experimentally  observed  by  White,  and  the  inverted  mode  was  again 
reported  to  give  lower  vibration  amplitudes  and  transraissibilities.  The 
performance  of  the  squeeze  film  was  shown  to  be  depedant  on  the  absolute 
values  of  a  particular  non-dimensional  bearing  parameter  and  on  a  shaft 
unbalance  parameter.  They  suggest  for  low  transmissibilities,  a  low 
centralising  spring  rate  should  be  used. 

Rabinowitz  and  Hahn  (24,  25,  26)  conducted  theoretical  and  experimental 
i lives igat ions  into  the  stability  of  a  centralised  squeeze  film  mounted 
flexible  rotor  with  rolling  element  healings.  They  show  that  the  squeeze 
film  pressurization  eliminates  bi-stable  operation,  the  lower  eccentricity 
mode  being  adopted,  and  allows  operation  through  the  first  pin  pin  critical 
speed  which  produced  violent  vibration  when  the  film  was  underpressurised, 
both  pressurization  (which  develops  a  2n  film)  and  an  increase  of  the 
bearing  parameter  discussed  in  reference  (23)  resulted  in  a  more  heavily 
damped  system.  line  tuning  could  be  achieved  by  controlling  the  lubricant 
temperature  (and  hence  its  viscosity),  the  control  of  which  was  found  not 
to  he  critical  for  close  to  optimum  performance.  Vibration  amplitude  at 
the  system  critical  speed  was  shown  to  be  attenuated  by  a  factor  of  up  to 
25  in  some  cases.  The  pressurization  levels  required  for  the  development 
of  a  2n  film  were  shown  to  be  as  had  been  earlier  predicted  by  Simandiri 
and  Hahn  (29). 
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In  1981  Hahn  (27)  developed  the  theoretical  model  described  in  reference 
(23)  to  allow  for  pressurization  resulting  in  a  full  277  film.  The  paper 
provides  design  data  for  effective  squeeze  film  performance  as  vibration 
isolators.  The  model  suggests  that  transmissibi lities  below  one  can  only 
be  achieved  in  certain  cases  at  the  expense  of  increasing  rotor  vibrational 
amplitudes.  The  bi-stable  operation  is  again  prevented  by  pressurization 
and  with  no  pressurization  there  is  a  critical  unbalance  above  which  the 
squeeze  film  effectiveness  is  minimal.  In  1976  Gunter,  Barrett  and  Allaire 
(31)  used  Reynolds  equation  to  produce  film  force  coefficients  for  each 
given  bearing  position  and  consequently  provide  data  for  optimum  damper 
design  in  respect  of  stability  and  transmissibi 1 ity.  Their  model  combines 
the  single  mass  rotor  and  bearing  equations  and  is  integrated  forward  in 
time  to  give  system  response  orbits.  Only  positive  pressures  were 
considered  to  contribute  to  film  forces.  They  showed  that  the  2tt  film 
reduces  the  stiffness  of  the  film  to  zero  but  doubles  the  film  damping  as 
compared  with  that  of  a  eavitated  (77)  film. 

Numerous  computet  predictions  of  shaft  trajector ies  were  plotted  and  the 
model  has  been  applied  to  several  compressors  with  successful  suppression 
of  non  syncronous  components  of  whirl.  They  conclude  that  with  short  rigid 
rotors,  reducing  squeeze  film  clearance  will  improve  stability  but.  for  long 
flexible  rotors  the  same  technique  can  be  catastrophic.  Increasing  the 
damping  and  effective  shaft  stiffness  may  also  improve  stability.  These 
conclusions  can  be  related  to  Gunter’s  work  of  1966  on  support  stiffness 
(13)  excessive  squeeze  film  stiffness  or  damping  were  both  shown  to  reduce 
damper  effectiveness  by  increasing  transmitted  forces.  Optimum  stability 
and  transmissibi 1 ity  values  were  for  eccentricity  ratios  below  O.d. 


Attention  was  later  focused  on  more  realistic  rotors  by  Cookson  and  Xin  Hai 
Feng  (32)  who  produced  an  assymetric  multi-mass  model  of  a  flexible  rotor 
supported  by  uncentralised  and  unpressurized  squeeze  films.  This  model 
gave  a  more  realistic  representation  of  a  practical  rotors  system  than  the 
usual  single  mass  model.  The  use  of  uncentralised  squeeze  films  means  that 
the  circular  syncronous  orbit  assumption  cannot  be  used  and  the  linear 
model  restricts  its  applicabilty  to  low  vibrational  amplitudes.  The 
analysis  shows  that  the  gravity  effects  increases  the  transmissibi 1  ity  of 
the  uncentralised  film  and  that  low  bearing  operating  parameter  values 
reduce  the  effect  of  foundation  stiffness  on  transmissibility,  amplitude 
ratio  and  system  critical  speed.  It  was  shown  that  the  transmitted  force 
could  be  kept  low  over  the  entire  speed  range,  and  by  varying  the  support 
stiffness  the  systems  critical  speed  could  be  shifted  in  the  frequency 
range. 

None  of  the  work  was  experimentally  verified  and  therefore  conclusions  that, 
have  not  been  previously  experimentally  confirmed  by  other  authors  must  be 
held  in  question. 

Theoretical  and  experimental  work  on  squeeze  films  without  a  centralising 
spring  for  both  rigid  and  flexible  rotors  was  published  by  Cookson  and 
Kossa  in  19B0  (33).  As  no  film  pressur i  sat i on  was  used  the  tt  film  case 
applies.  The  results  show  as  in  reference  (32)  that  the  gravity  parameter- 
should  be  kept  low  and  that  the  bearing  operating  parameter  discussed  in 
reference  (23)  is  optimum  at  0.1.  Instability  resulted  from  lower- 
operating  parameter  values  and  the  film  tended  to  behave  as  a  rigid  support 
at  higher  values.  The  authors  state  that  squeeze  films  do  not  reduce 
vibrational  effects  by  shifting  the  resonance  away  from  operating  speed  but 
that  they  reduce  the  amplitude  of  vibrations  by  a  vibration  redistribution 
process.  However  this  latter  conclusion  contradicts  the  theoretical  work 
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presented  in  reference  (32)  and  later  the  experimental  results  shown  in 
reference  (47). 

From  the  above  studies  the  ability  of  the  squeeze  film  to  attenuate 
vibration  amplitude  and  force  transmission  to  supports  was  shown  provided 
the  bearings  were  not  situated  at  rotor  free-free  nodes  and  that  the  film 
stiffness  was  not  such  as  to  promote  instability.  The  majority  of 
theoretical  models  showed  good  agreement  with  experimental  results  in  the 
lowe  eccentricity  rn age  but  tended  to  give  poor  prediction  accuracy  at  high 
(E>0.4)  eccentricity  ratios.  Work  by  Bannister  (15)  in  1972  on  the 
non-linear  characteristics  of  hydrodynamic  film  forces,  showed  the 
excessive  inaccuracies  of  linear  models  for  Ii>0.4.  By  developing  the 
second  differential  of  the  fluid  film  forces  with  respect  to  deflection, 
(ignored  by  linear  models)  he  produced  28  coefficients  to  model  the  film 
and  with  experimental  support,  showed  that  for  amplitudes  above  15%  of 
radial  clearance  greatly  improved  predictions  result. 

Despite  the  accepted  predominance  of  non-linear  oil  film  forces  at  high 
eccentricity  ratios,  the  simplification  of  the  analysis  with  the  consequent 
reduction  in  computer  time  required  still  leads  many  investigators  to  adopt 
linear  models.  If  the  orbit  amplitudes  are  kept  low  these  models  are 
sufficiently  accurate  to  justify  the  simplification. 


HYDROSTATIC  BEARINGS 

The  static  properties  of  hydrostatic  bearings  are  well  documented  by 
authors  Raimond  and  Boyd  (36)  and  with  additional  design  data  by  O’Donoghue 
and  Rowe  (35,46).  These  cover  investigations  into  various  betaring 
geometries  with  effects  on  static  loads  bearing  stiffnes,  oil  flow  rate 


etc,  usuully  the  restrictors  used  are  of  the  capillary  or  orifice  design. 


r 


In  1976  Singh  and  Singh  (37)  developed  a  method  of  stiffness  optimisation 
of  a  hydrostatic  bearing  by  the  use  of  spring  loaded  piston  type  variable 
restrictor.  This  was  similar  in  operating  principle  to  the  spring  loaded 
diagphram,  where  the  recess  back  pressure  forces  the  piston  back  resulting 
in  a  reduced  resistance  to  flow  and  an  increase  in  flow  rate  and  fluid 
pressure.  Unfortunately  Singh  gave  no  static  or  dynamic  performance  data 
for  his  design  but  did  provide  a  useful  demonstration  of  how  restrictor 
,  design  can  take  into  account  bearing  parameters  to  give  the  required 

|  stiffness. 

i 

The  papers  dealing  with  the  dynamic  characteristics  of  hydrostatic  bearings 
are  scarce  presumably  due  to  the,  to  date,  satisfactory  dynamic  performance 
of  hydrostatic  bearings  based  on  static  considerations.  However 
significant  contributions  have  been  made  by  the  following  authors. 

I.icht  and  Cooley  (38)  theoreticlaly  studied  the  dynamic  performance  of  u 
hydrostatic  bearing  and  concluded  that,  for  incompressible  flow,  large 
squeeze  numbers  (the  ratio  of  volume  flow  rate  due  to  vibration  to  static 
volume  flow  rate)  gives  high  damping  and  low  squeeze  numbers  gives  low 
damping.  Unfortunately  the  effects  of  restrictor  type  or  vibrational 
frequency  are  not  considered. 

lngert  (39)  in  1973  developed  a  bearing  characteristic  equation  and  used 
its  roots  to  define  the  limits  of  stability.  He  concludes  that  for 
realistic  pad  geometries  instability  is  unlikely  and  that  high  oil 
viscosities,  high  bearing  area  and  and  low  bearing  clearance,  and  for 
compressible  flow,  higher  pocket  pressure  results  in  reduced  vibrational 
amplitude.  Each  of  these  factors  leads  to  higher  bearing  stiffness. 
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It  was  shown  by  Newton  and  Howarth  in  1974  (4)  that  the  tilt  and  dynamic 
action  of  a  hydrostatic  bearing  would  change  the  static  operating  variables 
and  conclude  that  static  stiffness  and  damping  are  highest  for  no  tilt. 
Unfortunately  the  model  does  not  consider  dynamic  stiffness  and  damping  and 
ignores  the  effect  of  inertia  by  the  use  of  steady  velocities.  Cooms  and 
Dowson  in  1965  (47)  published  work  showing  that  for  static  considerations 
inertia  effects  were  insignificant  except  at  very  low  suppiy  pressures. 
The  importance  of  this  work  is  that  it  provided  qualitive  data  to  support 
the  assumptions  made  in  most  investigations  of  the  static  properties  of 
hydrostatic  bearings,  that  the  inertia  effects  can  be  neglected.  However, 
the  effect  of  inertia  on  the  dynamic  properties  of  u  hydrostatic  bearing 
with  a  relatively  large  restrictor  fluid  volume  can  still  be  considered 
significant  (49). 

Brown  in  1961  (41)  analyses  the  dynamic  behaviour  of  hydrostatic  bearings 
by  assuming  compressible  laminar  flow  and  neglecting  fluid  inertia.  He 
represents  the  bearing  system  by  a  series  damper  and  spring  to  represent 
the  dynamic  stiffness  and  damping  with  a  parallel  spring  to  represent  the 
static  stiffness,  hence  assuming  that  they  are  independent.  He  suggests 
that  the  dynamic  properties  were  largely  dependent  upon  the  fluid  under  the 
sills,  i.e.  the  fluid  in  the  pocket  and  between  the  pocket  and  restrictor 
could  be  neglected  when  deriving  compressibility  and  that  the  damping  was 
independent  of  inlet  restrictor.  These  assumptions,  although  shown  to  be 
invalid  (48)  gave  good  results  up  to  20  Hz  but  large  errors  above  this 
frequency.  He  concludes  that  dynamic  stiffness  is  always  higher  than 
static  stiffness  and  that  damping  is  always  high. 


13- 


HYDROSTATIC  SQUEEZE  FILM 


Due  to  the  need  for  stability  of  hydrostatic  journal  bearings  of  higher 
opearating  speed,  Choy  and  Halloran  in  1982  (42)  conducted  a  theoretical 
and  experimental  investigation  into  the  damping  and  stabilising  effects  of 
a  hydrostatic  journal  bearing.  The  linear  damping  forces  are  calculated 
using  Reynolds  equation  and  the  equilibrium  positon  found  by  iteration  of 
pocket  pressure  by  equating  the  inlet  and  leakage  flow  of  each  pocket.  The 
film  showed  very  little  change  in  stiffness  or  damping  with  supply  pressure 
increase  but  both  the  stiffness  and  damping  of  the  film  changed 
significantly  with  eccentricity  variation;  the  maximum  values  occuring 
higher  eccentricity  levels  for  larger  film  clearances  but  the  peak  values 
reducing  in  amplitude  as  clearance  increased.  The  damping  is  shown  to  be 
due  to  the  fact  that  the  oil  must  pass  over  the  sill  and  is  not  dependent 
on  the  amount  of  flow  to  the  recess. 

Koshal  and  Rowe  in  1981  (44,45)  experimentally  and  theoretically 
investigated  the  hybrid  operation  of  hydrostatic  bearings  and  confirmed  the 
benefits  of  high  load  capacity  at  zero  speed  and  high  stiffness  for  low 
eccentricity  operation,  attributable  to  hydrostatic  action  and  high 
stiffness  and  load  capacity  at  high  eccentricity,  attributable  to 
hydrodynumic  action.  The  experimental  runs  shown  a  possibility  (at  high 
power  rates,  frictional  power/pumping  power,  i.e.  E>0.75)  of  negative 
bearing  stiffness  (oil  film  collapse).  They  suggest  that  this  is  due  to 
high  film  thickness  gradients  at  high  eccentricity  ratios  with  the 
associated  circumferential  flow  and  severe  circumferential  pressure 
gradients  with  the  resultant  reduced  effective  high  pressure  area.  The 
increase  in  slot  pressure  as  eccentricity  tends  to  one  cannot  compensate 
for  the  loss  in  effective  high  pressure  area  and  the  bearing  collapses. 
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As  the  phenomena  was  never  experienced  during  their  experiments,  it  is  hard 
to  know  whether  the  explanation  is  correct  or  indeed  whether  the  phenomena 
exists. 

VARIABLE  IMPEDANCE  BEARINGS 


The  dynamic  behaviour  of  a  shaft  bearing  system,  as  with  ail  mechanical 
systems,  is  determined  by  the  mechanical  impedances  or  more  accurately  the 
interaction  of  the  mechanical  impedances  of  the  system.  These  impedances 
can  be  represented  in  complex  form  with  the  intertias  and  stiffnesses 
represented  by  opposing  imaginary  components  and  the  damping  represented  by 
the  real  component.  The  impedances  are  a  function  of  frequency  and  at 
certain  frequencies  the  imaginary  components  cancel  out  and  system 
instability  results. 

The  above  suggests  that  the  dynamic  behaviour  of  a  rotating  shaft  system 
can  be  controlled  by  the  control  of  the  system  impedances.  Due  to  the 
limitations  put  on  turbine  designers,  the  inass,  internal  damping  and 
stiffness  of  the  components  are  generally  fixed  for  a  particular  system, 
the  support  structure  being  the  only  component  available  for  modification 
of  impedance  values. 

Goodwin  in  1981  (43)  experimentally  and  theoretically  investigated  a 

hydrodynamic  journal  bearing  surrounded  by  an  outer  hydrostatic  bearing 
which  was  pinned  to  prevent  rotation.  The  hydrodynamic  film  was  modelled 
using  eight  linear  film  coefficients  and  the  effect  of  variation  of  the 
bearing  parameters  on  the  system  performance  was  investigated.  The  theory 
and  experimental  results  show  the  ability  to  independently  control  the 
hydrostatic  static  stiffness  by  altering  the  setting  of  a  flow  restrictor 
in  the  hydrostatic  film  supply  line,  whilst  the  X  and  Y  axis  dynamic 
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stiffness  and  damping  can  be  controlled  by  varying  the  flow  resistance  of  a 
restrictor  which  linked  the  hydrostatic  bearing  pockets  to  an  accumulator. 
Substantial  variations  of  the  overall  system  stiffness  and  damping  were 
recorded  due  to  the  introduction  of  hydrostatic  film.  These  were  measured 
at  various  excitation  frequencies.  Of  particular  reference  was  the  shaft 
of  system  critical  speed  by  350  rpm  and  the  reduction  of  vibrational 
amplitude  by  80  per  cent  predicted  for  a  large  machine  when  operated  with 
hydrostatic  supports. 

Burrows,  Sahinkaya  and  Turkay  (49)  in  1984  carried  out  theoretical  and 
experimental  investigations  into  the  effect  of  optimising  the  damping  for  a 
given  support  stiffness  and  shaft  speed.  •  The  damping  was  introduced  by  a 
squeeze  film  surrounding  a  rolling  element  bearing  and  controlled  by 
applying  a  step  change  in  the  oil  supply  pressure.  As  shown  by  (42)  no 
damping  changes  occur  with  supply  pressure  increase  unless  the  pressure 
increase  is  supressing  the  onset  of  cavitation  in  the  squeeze  film. 
Presumably  this  control  of  cavitation  is  the  mechanism  being  utilised  in 
this  investigation.  They  state  in  conclusion  that  for  the  conditions  of 
their  experiments  the  adaptive  dampers  had  little  effect  on  amplitude 
compared  with  passive  dampers  when  the  squeeze  f i 1ms  were  situated  at  the? 
shaft  ends  and  that  the  transmitted  force  increased. 

They  do  however  display  the  beneficial  effect  of  Belleville  washers  in 
changing  the  support  stiffness  and  being  able  to  ’tune’  the  system  support 
damping  for  a  given  support  stiffness  which  together  modify  the  critical 
speed  and  response  amplitude. 
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In  1979  Sandler  (48)  investigated  the  automatic  control  of  the  dynamic 
response  of  a  rotating  shaft  system.  The  experiments  were  based  on  the 
variation  of  the  shaft  stiffness  by  automatically  controlling  the  bearing 
position  on  the  shaft.  Two  bearing  positions  were  possible  each  with  an 
associated  resonance  peak  for  each  critical  speed.  The  result  of  changing 
the  bearing  position  during  speed  increase  or  decrease  was  that  the  system 
adopted  the  response  of  the  associated  bearing  position  i.e.  the  lower  of 
the  two  dynamic  responses.  Despite  the  obvious  limitations  of  changing 
shaft  stiffness  on  a  practical  system,  the  work  does  show  how  automatic 
control  of  system  impedances,  in  this  case  shaft  stiffness,  can  be  utilised 
to  modify  the  dynamic  behaviour  of  a  system. 

The  increase  in  power  output  of  turbine  systems  over  the  past  thirty  years 
has  been  accompanied  by  higher  running  speeds  and  rotor  length  to  diameter 
ratios,  with  the  result  that  systems  often  operate  above  several  of  their 
critical  speeds.  Because  of  the  difficulty  in  accurately  predicting 
support  dynamic  impedances  and  rotor  dynamic  behaviour,  rotor  system 
designers  sometimes  find  predictions  of  rotor  system  critical  frequencies 
sufficiently  inaccurate  as  to  result  in  the  operating  speed  coinciding  with 
or  close  to  the  system  critical  speed.  This  has  resulted  in  high 
vibrational  amplitudes  and  transmitted  forces  at.  the  normal  operating 
speed. 

The  ability  to  control  the  support  impedunec  of  the  shaft  system  after  the 
system  has  been  installed  will  result  in  a  control  of  system  critical 
frequencies  with  the  consequent  ability  to  shift  them  away  from  operating 
speeds  and  the  resultant  reduction  in  response  amplitudes.  For  system 
speed  run  up  or  for  a  variable  speed  system  the  ability  to  control  the 
support  impedance  (particularly  stiffness)  over  a  range  of  values  depending 
on  shaft  speed,  makes  available  the  possibility  of  a  speed  controlled 
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automatic  system  which  will  produce  low  response  amplitudes  at  any  running 
speed . 

The  possibility  of  variable  impedance  support  control  system  failure  must 
also  be  considered  during  the  development  of  the  system. 
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STATIC  CHARACTERISTICS 


SINGLE  HYDROSTATIC  PAD  AND  FULL  HYDROSTATIC  RING 
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STATIC  CHARACTERISTICS  OF  A  SINGLE  CURVED  HYDROSTATIC  PAD 
INTRODUCTION 

The  proposed  solution  to  the  requirement  for  a  variable  impedance  support 
for  turbine  systems  is  based  on  a  modified  hydrostatic  bearing  system.  The 
main  bearing,  required  to  accoraadate  shaft  rotation,  will  be  a  rolling 
element  type,  with  the  outer  surface  being  prevented  from  rotating  and 
acting  as  one  surface  of  the  hydrostatic  pads.  The  bearing  outer  surface 
will  however  be  free  to  move  within  the  hydrostatic  support  clearance.  As  a 
result  of  this  relative  motion,  together  with  a  continuous  oil  supply  to 
the  clearance,  the  land  areas  of  the  hydrostatic  support  will  act  a  partial 
squeeze  film,  sometimes  used  in  its  full  form  to  introduce  system  damping. 
The  hydrostatic  effect  will  supply  the  static  stiffness  to  oppose  the 
gravity  load  of  the  shaft  system  and  centralise  the  shaft,  this  is  normally 
supplied  by  a  form  of  mechanical  centralising  spring  with  the  full  squeeze 
film  systems.  As  a  result  of  this  inherent  centralising  stiffness  from  the 
hydrostatic  effect,  the  design  will  possess  an  oil  filled  recess.  The 
presence  of  this  recess  results  in  reduced  land  area  and  hence  reduced 
sqoeze  film  forces  but  also  introduces  the  possibility  of  varying  the 
dynamic  stiffness  of  the  support  independant ly  to  the  static  stiffness. 

The  dynamic  characteristics  of  a  hydrostatic  system  are  dominated  by  the 
effect  of  flow  variations  imposed  by  the  dynamic  load.  When  the  load 
approaches  the  hydrostatic  pad,  the  recess  pressure  will  increase  and  the 
supply  flow  reduce  and  in  certain  cases  reverse.  The  flow  variations  can  be 
controlled  by  incorporating  an  additional  exit  for  the  recess  oil.  An 
accumulator  can  therefore  be  introduced  into  the  system,  the  effect  of 
which  can  be  varied  by  varying  the  value  of  the  assosiated  accumulator 
restrictor.  With  this  modification,  as  the  main  bearing  outer  shell 
approaches  the  hydrostatic  recess,  the  recess  oil  will  be  partially 
compressed  and  partially  accelerated  to  flow  either  over  the  lands,  down 
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the  supply  restrictor  or  down  the  accumulator  restrictor.  It  will  therefore 
be  possible  to  control  this  oil  flow  and  consequently  modify  the  supports 
displacement  opposing  forces  (stiffness)  flue  to  recess  fluid 
compressibility  and  the  velocity  opposing  forces  (damping)  thus  meeting  the 
requirement  for  a  variable  impedance  support  to  modify  the  rotor  systems 
dynamic  response. 


SINGLE  PAD  ANALYSIS 

The  analysis  that  follows  is  in  reference  to  a  single  curved  hydrostatic 
pad.  The  responses  given  are  those  due  to  vertical  deflection  of  the 
supported  load  with  velocity  and  acceleration  both  zero.  The  dimensions 
and  other  pad  values  are  similar  to  those  of  practical  use  for  a  1/20  model 
of  a  typical  750  kN  rotor  system  supported  with  one  pad  at  each  end,  the 
theory  is  however  applicable  to  any  pad  dimensions  with  any  exceptions 
commented  on  in  the  text. 

To  allow  simplification  of  the  analysis,  the  leakage  pressure  of  the  pad, 
i.e.  the  fluid  pressure  at  the  land  outer  perimeter,  is  assumed  zero 
(relative  to  atmospheric  pressure).  When  the  multi  pad  system  is  analysed 
this  enables  inter  pad  flow  and  other  pad  interactions  to  be  neglected. 
The  practical  realisation  of  no  inter  pad  flow  is  achieved  by  drainage 
channels  running  along  the  axial  length  of  the  hydrostatic  ring  between 
each  of  the  pads. 

The  system  can  be  considered  as  one  fixed  and  one  variable  resistance 
representing  the  input  restrictor  and  the  parallel  plate  land  resistance 
respectively.  A  static  supply  pressure  causes  a  flow  through  tl  .  pad  the 
magnitude  of  which  is  controlled  by  the  two  resistance  values.  The 
pressure  drop  ucross  the  input  restrictor  is  directly  proportional  to  the 
volumetric  flow  rate  and  thus  determines  the  pocket  pressure  which  opposes 

21  . 


the  load.  The  parallel  plate  resistance  is  proportional  to  the  inverse  of 
the  land  clearance  cubed  which  in  series  with  the  restrictor  resistance 
determines  the  flow  rate.  Consequently  by  conservation  of  flow  the  pocket 
pressure  can  be  determined  for  a  given  pad  deflection. 

The  electrical  analogy  of  the  system  is  a  supply  voltage,  a  fixed 
resistance,  a  load  voltage  and  a  series  resistance  a  1/  charge^,  the 
current  representing  the  volume  flow  rate.  Using  this  it  can  be  seen  Fig. 
1  that  the  resistances  are  acting  as  dampers  in  the  system  i.e.  generating 
a  back  pressure  proportional  to  the  flow  rate.  Despite  this  damping  system 
the  overall  effect  is  to  change  the  pocket  pressure  which  will  be  ’seen’  as 
a  stiffness  by  the  applied  load.  This  apparent  90°  phase  shift  of  the 
fluid  flow  damping  to  generate  a  pad  stiffness  effect  (Fig.l)  is  due  to  the 
steady  fluid  flow  required  to  support  the  load.  Thus  this  fluid  flow 
damping  can  be  substituted  by  a  spring  in  the  mechanical  analogy  whose 
stiffness  is  constant  (Linear  model)  or  f(e)  non-linear  model.  It  should 
be  noted  that  the  fluid  compressibility,  system  resistance  and  fluid 
inertia  do  not  undergo  similar  90°  phase  shifts  for  the  dynamic  model  and 
remain  stiffness,  damping  and  inertia  respectively. 
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MIDLAND  APPROXIMATION 


The  land  parallel  plate  resistance  causes  a  linear  pressure  drop  with 
distance  across  the  land,  as  the  pressure  drop  is  linear  then  this  pressure 
gradient  can  be  substituted  for  a  step  drop  in  pressure  at  the  end  of  a 
land  of  half  the  original  width.  Thus  this  midland  approximation  can  be 
used  to  simplify  the  complex  pressure  profile  of  the  hydrostatic  pad. 
Goodwin  (42)  shows  that  for  a  pad  recess  pressures  of  10  Pa  and  a 
consequent  pad  load  of  17.7  x  10"^  N  an  error  of  1.7%  above  true  load  is 
calculated  using  this  midland  approximation.  Due  to  the  more  rapid 
pressure  drop  with  distance  at  the  land  corners  the  midland  accuracy  will 
be  dependent  on  the  percentage  of  land  corner  area  to  straight  land  area, 
the  example  given  by  Goodwin  being  for  a  percentage  corner  to  total  land 
are  of  22%.  Any  significant  increase  above  this  value  may  require  a 
correction  factor  for  the  final  model. 

DISCUSSION 

The  graph  of  load  v  eccentricity  (appendix  1)  shows  the  trend  for  pad  load 
to  tend  to  a  maximum  value  of  (Ps  xAe)  as  eccentricity  reduces  to  "1.0  for 
all  B  values.  This  trend  can  be  explained  by  referring  to  the  mechanics  of 
the  pad  system. 

As  eccentricity  approaches  "1.0  the  land  clearance  tends  to  zero  and  the 
associated  resistance  tends  to  infinity.  The  resultant  zero  flow  rate 
causes  the  pressure  drop  across  the  inflow  restrictor  to  become  zero  and 
the  pad  pressure  ,  tends  to  the  supply  pressure  independent  of  inflow 
resistance  and  therefore  P  value. 
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As  the  eccentricity  increases  above  £  =  “1.0  the  land  resistance  assumes  a 
finite  value  and  flow  recommences.  The  rate  at  which  the  pressure  drop 
across  the  restrictor  increases  with  flow  rate  increase  is  by  definition 
the  restrictor  resistance  value,  thus  the  higher  the  restrictor  value  the 
faster  the  pressure  drop  increases  with  flow.  This  is  shown  by  graph  1 
with  the  lower  P  values  displaying  higher  load  eccentricity  gradients. 

The  stiffness  of  the  pad  is  defined  by  the  rate  of  change  of  load  with 
deflection  of  the  shaft  which  is  given  by  the  gradient  of  the  load 
eccentricity  graphs.  As  shown  by  graph  2,  the  magnitude  of  maximum 

stiffness  and  the  eccentricity  at  which  it  occurs  increase  with  P  ratio 
reduction,  these  maximum  stiffness  values  represent  the  points  of 
inflection  of  the  load  curves.  The  lowest  stiffness  values  occur  at  an 
eccentricity  value  of  e  =  “1.0  due  to  the  asymptotic  approach  of  the  load 
curves  to  the  limiting  value.  Work  by  Burrows,  Sahinkaya  and  Turkay  (49) 
and  Goodwin  (42)  show  that  the  first  and  second  critical  speeds  of  a 
typical  system  are  most  sensitive  to  support  stiffness  changes  when  the 
support  stiffness  is  of  the  order  of  5  x  10'^  NM  ^  and  that  the  third  and 
fourth  critical  speeds  are  most  sensitive  to  support  stiffness  changes  when 
the  support  stiffness  is  of  the  order  of  1  x  10^  NM“^.  Thus  a  range  from 
5  x  1()4  to  1  x  10^  NM“1  would  be  desirable  in  shifting  the  critical  speeds 
of  the  system  away  from  running  speed. 

To  achieve  this  low  value,  the  stiffness  graph  suggests  operating  at  high  P 
values  (approaching  1.0)  and  eccentricities  below  £  =  “0.9  or  lower  P 
values  (below  0.1)  and  eccentricities  above  £  =  0.9. 
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Adopting  the  high  P  solution  would  require  high  static  deflections  which 
may  impose  unacceptably  high  stress  levels  on  the  rotor  system  due  to  shaft 
misalignment.  Uncentralised  squeeze  films  do  however  successfully  operate 
at  similarly  high  eccentricity  ratios  particularly  at  the  low  speed  range 
which  suggests  that  the  system  could  accommodate  such  high  eccentricity 
values.  However,  referring  to  graph  1,  this  shows  that  this  would  require 
a  reduction  in  supply  pressure  to  give  the  required  working  load  which 
would  result  in  no  overload  capability.  Using  J3  =  U. 1  necessitates 
eccentricity  values  of  €  =  >+0.9.  Referring  to  graph  1  the  supply  pressure 
increase  required  would  result  in  unacceptably  high  flow  rates.  The  pad 
stiffness  is  also  proportional  to  supply  pressure  and  consequently  an 
increase  in  supply  pressure  would  also  result  in  a  corresponding  stiffness 
increase  thus  the  low  P  solution  is  also  unacceptable. 

Referring  to  equation  (7)  the  final  possibilities  is  an  increase  of  the 
land  clearance  at  e  =  0.  This  however  both  increases  the  flow  rate  (graph 
4)  and  the  deflection  required  before  load  increases  are  opposed  by  the 
pad,  consequently  shaft  stresses  would  increase.  Increase  of  the  clearance 
does  however  reduce  the  pad  stiffness  by  reducing  the  rate  of  flow  change 
w.r.t.  deflection.  Unfortunately  (see  gruph  4)  the  associated  flow 
increase  is  severe  enough  to  make  clearance  increase  unsuitable  as  a  means 
of  reducing  pad  stiffness. 

The  above  analysis  shows  that  modification  of  standard  hydrostatic  pad 
parameters  is  not  a  suitable  method  of  reducing  support  stiffness  to  the 
value  required  and  that  the  method  of  accumulator  installation  must  be 
investigated  in  order  to  change  the  dynamic  stiffness  independently  of  the 
static  stiffness. 
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The  full  hydrostatic  support  consists  of  four  pads  at  90  degree  intervals 
around  the  shaft  axis  and  in  the  same  plane  as  the  main  shaft  bearings. 
Drainage  channels  at  atmospheric  pressure  run  axially  between  each  pad  and 
prevent  pad  flow  interaction.  The  characteristics  of  each  pad  can  therefore 
be  calculated  individually,  with  consideration  given  to  the  different 
relative  shaft  position,  and  the  results  directly  added  together. 

As  the  pads  are  identical  and  the  pad  system  symetrical  about  both  the  X 
and  Y  axis,  the  graphs  for  the  four  pad  system  show  the  same  symetry 
(Graphs  6-11).  The  pad  reaction  force  decreases  with  shaft  distance  from 
the  pad,  however  the  opposing  pad  force  will  increase  until  at  the  shaft 
central  position,  the  two  reactions  will  cancel  giving  a  zero  resultant 
force.  As  a  result  of  this,  despite  the  variations  in  pad  normal  reaction 
force  as  the  shaft  moves  across  the  pad  face,  the  resultant  force  remains 
zero  (Graph  10).  This  is  not  the  case  for  the  pad  system  stiffness  profile 
(Graph  9).  The  static  stiffnesses  of  the  pads  arc  always  positive  and  the 
summation  of  such  stiffnesses  always  give*s  a  positive  resultant.  For  fl>, 0.5 
the  maximum  stiffness  uecurs  at  the  central  shaft  position  (Graph  7)  due 
to  the  minimal  peak  and  gradual  stiffness  changes  with  normal  shaft 
deflections  (Graph  2).  For  lower  P  values,  the  maximum  stiffness  peak  for 
the  single  pad  (Graph  2)  is  severe  enough  to  remain  when  the  opposing  pad 
effects  are  added. 

Despite  the  lack  of  reaction  force  from  the  land  surfaces,  no  accut.o 
reduction  in  reaction  force  accurs  when  the  shaft  is  loaded  in  a  direction 
along  the  shaft  centre  to  drainage  channel  centre  line  with  the  pad 
dimensions  used  (Graph  13).  This  suggests  that  no  extreme  radial  static 
stiffness  variations  will  occur  as  the  shaft  travels  around  the  hydrostatic 
clearance  circle  due  to  shaft  unbalance,  even  if  eccentricities  are  close 
to  unity.  This  situation  must  be  reassessed  for  the  final  pad  dimensions  to 
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ensure  that  minimal  angular  variations  still  exist. 

Graph  12  shows  the  cross  coupled  stiffnesses  for  the  four  pad  system  ie. 
the  effect  on  stiffness  characteristics  in  one  direction  due  to  shaft 
motion  in  the  normal  direction  in  the  bearing  plane.  As  shown,  the  cross 
coupled  effects  remain  zero  when  the  shaft  is  on  axis  X  or  Y.  This 
simplifies  the  situation  of  a  static  load  acting  through  the  lower  pad 
upper  pad  centre  line.  However,  the  off  axis  cross  couplings  are  non  zero 
and  must  therefore  be  taken  into  account  if  either  the  static  loading  is 
not  vertical  or  if  the  position  of  the  shaft  changes  in  the  bearing  plane 
due  to  unbalance  forces  lifting  and  orbiting  the  main  bearings. 

The  flow  rates  of  the  total  support  are  simply  the  sum  of  the  individual 
pads  calculated  with  the  shaft  in  its  correct  poition  relative  to  the  pad. 
The  shaft  centred  flow  rate  is  therefore  four  times  the  single  pad  e=0  flow 
rate  (Graph  3  &  11).  For  P=0.9  as  the  shaft  shall  moves  from  e=0  away  from 
the  pad,  the  flow  rate  increases  far  more  sharply  than  the  flow  from  the 
oppisit  pad  decreases,  consequently  the  flow  rate  increases  with 
eccentricity  ratio  for  high  P  values.  For  /3=0.i  the  flow  variation  with 
eccentricity  is  very  small  and  consequently  the  sum  of  the  opposing  pad 
flows  remains  almost  constant.  Obviously  the  magnetude  of  the  high  P  value 
flows  are  higher  due  to  the  higher  pocket  pressure  for  a  given 
eecentrici ty.  Due  to  the  predictability  of  the  flow  patterns  no  problems 
should  occur  as  a  result  of  excessive  support  flow  rates. 

The  above  analysis  is  not  altered  by  the  additon  of  an  accumulator  and  the 
assosiated  resistance  Raccum  (Fig.l)  regardless  of  the  accumulator 
operating  parameter  AP/  AVol.  The  flow  into  the  accumulator  will  stop 
when  the  accumulator  pressure  equals  the  pocket  pressure.  Thus  the  static 
stiffness,  load  and  flow  remain  unaltered.  The  R  value  can  be  assumed  zero 


for  small  amplitude  dynamic  variations  but  obviously  has  a  finite  .  It  can 
therefore  be  represented  by  a  capacitor  in  the  electrical  network  Fig  1. 
With  a  change  in  load  voltage  this  capacitor  charges  but  the  network  is 
subsequently  unaffected. 
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The  above  was  converted  into  computer  program  and  the  results  plotted. 
Appendix. 
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